function [ig,ie,xi] = locate_point(x,grid,ddc_grid,varargin)
% [ig,ie,xi] = locate_point(x,grid,ddc_grid)
% [ig,ie,xi] = locate_point(x,grid,ddc_grid,ie)
% [ig,ie,xi] = locate_point(x,grid,ddc_grid,ie,dont_stop)
%
% Locate the point x in a given grid. The optional input arguments are
% ie: initial guess
% dont_stop: for missing points, do not print an error but return a
%            negative value
%
% It is possible to pass ddc_grid = [].

 toll = 1e-10;

 d = grid{1}.d;

 checked = zeros(max([[grid{:}].ne]),length(grid));
 for ig=1:size(checked,2)
   checked(grid{ig}.ne+1:end,ig) = 1;
 end

 if(not(isempty(ddc_grid)))
   ddc_marker = ddc_grid{1}.ddc_marker;
 else % pass ddc_marker = [] if not used
   ddc_marker = -1;
 end

 found = 0;
 if(nargin>3)
   ig = varargin{1}(1);
   ie = varargin{1}(2);
 else
   ie = locate_candidate(checked);
 end

 while(found==0)

   % avoid infinite loops
   if(ie<0) % no candidates
     if(nargin>4)
       xi(1:d+1) = NaN;
       return
     else
       error(["Unable to locate point (",num2str(x'),")."]);
     end
   end
   e = grid{ig}.e(ie);

   xi(2:d+1) = e.bi*(x-e.x0); % barycentric coords
   xi(1) = 1 - sum(xi(2:d+1));
   [ximin iximin] = min(xi);
   checked(ie,ig) = 1;
   if(ximin>=-toll)
     found = 1;
   else
     ie = e.ie(iximin);
     found_new_ie = 0;
     if(ie>0)
       if(~checked(ie,ig))
         found_new_ie = 1;
       end
     else
       if(ie==-ddc_marker)
         ins = ddc_grid{ig}.gs(2,e.is(iximin));
         new_ig = ddc_grid{ig}.ns(2,ins)+1;
         ie = grid{new_ig}.s(ddc_grid{ig}.ns(3,ins)).ie(1);
         ig = new_ig;
         if(~checked(ie,ig))
           found_new_ie = 1;
         end
       end
     end
     if(~found_new_ie)
       [ig,ie] = locate_candidate(checked);
     end
   end
 end

return

function [ig,ie] = locate_candidate(checked)
 [ie,ig] = find(checked~=1);
 if(isempty(ie))
   ie = -1;
   ig = -1;
 else
   ie = ie(1);
   ig = ig(1);
 endif
return

